#! /usr/bin/perl -w 

######################################################################
# Time-stamp: <04/08/19 10:00:00 phlee>
# used for Package: EP::EPC::Core::Threeway_Similarity_Analysis
# Interpro Similarity Precomputation 
######################################################################


package main;

use strict;

use lib '/usr/local/EPNG/cgi-lib/perl/';
use EP::Config;


# constants
use constant MY_DEBUG => 0;
use constant MAIN_DIR => $EP::Config::h->{TSA_DATA_PATH};

# constants for &compute_batch_interpro_similarity(RESULT_DIR);
use constant RESULT_DIR => MAIN_DIR . "/interpro/result/";
use constant DATA_DIR => MAIN_DIR . "/interpro/data/";

use constant IDA_TERM_FILE_NAME => DATA_DIR . "IDA_term.txt";
use constant PROTEIN_IDA_FILE_NAME => DATA_DIR . "ProteinAC_IDA.txt";
use constant VECTOR_FILE_NAME => DATA_DIR . "vector.txt";

use constant PROTEIN2YEAST_FILE_NAME => MAIN_DIR . "/ac2gn/ac2gn.yeast";
use constant INTERPRO_FILE_EXT => ".summary";

# constants for &generate_IDA_img_URL();
use constant MOTIF_IMG_URL_FILE_NAME => MAIN_DIR . "/interpro/Protein_Motif_URL.txt";

# global variables 
my (@interpro_vector, %interpro_hash);
my (@protein_vector, @protein_interpro_array, $protein_interpro_2d_array);
my ($protein_cnt, $interpro_cnt);
my (%protein2orf);
my (%match_ratio_result);



#---------------------------------------------------------
# 1. INITIALIZATION
#    this should be called before using any one of functions
#    compute_batch_interpro_similarity();  or
#    generating IDA_img_URL();
#---------------------------------------------------------
    &init();

#-------------------------------------------------------------
# 2. COMPUTE THE SIMILARITY
#    for computing motif similarity between all pairs of genes
#-------------------------------------------------------------
    &compute_batch_interpro_similarity(RESULT_DIR);

#---------------------------------------------------------------
# 2. GENERATE the MOTIF IMAGE URL INFORMATION FILE
#    for loading motif image from INTERPRO website, 
#    whole url information is precomputed and stored in the file
#---------------------------------------------------------------
    &generate_IDA_img_URL(MOTIF_IMG_URL_FILE_NAME);



###################################################################################
# TimeStamp 2004-08-12 Phil Hyoun Lee
#
# METHOD     generate_IDA_img_URL 
# PURPOSE    for generating the motif url file 
# ARGS       $outputfilepath 
###################################################################################
sub generate_IDA_img_URL {

  print "\n\n********************************************\n";
  print "generate_IDA_img_URL() ....\n";
  print "********************************************\n";

  # open the output file
  my $filename = shift; 
  open FH, ">$filename";

  # for every protein,
  for(my $i=0; $i<$protein_cnt; $i++) {
        my $protein = $protein_vector[$i];  
        my $line = "http://www.ebi.ac.uk/interpro/matchline?code=A$protein,C,";

        # add distinct IPRO
        my $added_no=0;
        my $color_code=0;
        for(my $j=0; $j<$interpro_cnt; $j++) {
            if( $protein_interpro_2d_array->[$i][$j]==1) {
                if( $added_no > 0 ) { 
		    $line .= ",";
                }
                $line .= "IPR";

                # add zeros 
                my $interpro = $interpro_vector[$j];
                my $len = length($interpro);
                while ($len<6) {
			$line .= "0";
                        $len++;
                }  

		$line .= "$interpro:$color_code";

                $added_no++;
                $color_code++;
	    }
	}

	if( exists($protein2orf{$protein}) && $added_no>0 ) {
	    print FH "$protein2orf{$protein}\t$line\t$added_no\t$protein\n";
        }
  }

  close FH;
}

###################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     load_protein2orf 
# PURPOSE    for mapping the proteinAC to Yeast ORF, load the data   
# ARGS       $filepath 
# RETURN     ex. result is loaded into the global var. %protein2orf 
###################################################################################
sub load_protein2orf {

	my $filepath = shift;

	my (@list, $line, @protein_list, @orf_list, $protein, $orf);

	# open the file
	open FH, "<$filepath"
		or die "[ERROR] can't open the file $filepath : $!";

	# read the content
        while( defined($line=<FH>) ) {
                chomp($line);

                # if this is the comment or blank line, skip it.
		@list = split /\t/, $line;

		# check if protein has the correspondencing orfs.
		if( $#list==1 ) {
			@protein_list = split /,/, $list[0];
			@orf_list = split /,/, $list[1];

			# register it
			foreach $protein (@protein_list) {
				foreach $orf (@orf_list) {
					$protein2orf{$protein} = $orf;
				}
			}
		}

	}

	close FH;

	if(MY_DEBUG) { 
		while( ($protein, $orf) = each %protein2orf ) {
			print $protein, "=>", $orf, "\n";
		}
	}
} 


###################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     compute_batch_interpro_similarity 
# PURPOSE    compute_batch_interpro_similarity
# ARGS       $output_path;    ex. only the main directory
# RETURN     
###################################################################################
sub compute_batch_interpro_similarity {

        print "\n\n********************************************\n";
        print "compute_batch_interpro_similarity() ....\n";
        print "********************************************\n";

   	my ($output_dir, $output_path) = shift;	
	my $match_ratio;

	print "protein no $protein_cnt vector no $interpro_cnt\n";

	for(my $i=0; $i<$protein_cnt; $i++) {
		# this $protein_vector[$i] is the base protein

		# store the result BASE_PROTEIN.interpro_match.summary
		# but we should think about the ID issue here !
		# is it YEAST ORFs or PROTEIN ACs ?
		if( exists($protein2orf{$protein_vector[$i]}) ) {
			# if it has the matched ORF, store it with that name
			$output_path = $output_dir . $protein2orf{$protein_vector[$i]} . INTERPRO_FILE_EXT;
		}
		else {
                        # currently just skip it.
			next;

			# stored with protein name
			$output_path = $output_dir . $protein_vector[$i] . INTERPRO_FILE_EXT;
		}

		# match ratio is computed for all other proteins
		&init_match_ratio_result();
		for(my $j=0; $j<$protein_cnt; $j++) {

			$match_ratio = &compute_interpro_similarity($i, $j);	

			$match_ratio_result{$protein_vector[$j]} = $match_ratio;
		}

		# sort the tmp_result by the  match_ratio 
       		my @tmp_result = sort sort_by_match_ratio keys %match_ratio_result;


		# store the sorted result
		open FH, ">$output_path"
			or die "[ERROR] can't open the file $output_path : $!";
			

		my ($key, $value);
		print "$output_path processing ...........................\n";
               	for(my $tmp_index=$#tmp_result; $tmp_index>0; $tmp_index--) {
			$key = $tmp_result[$tmp_index];

			# if ORF exists
			if ( $match_ratio_result{$key} == -1 )  {
				# skip this 
				print "$key\n";
			}
		        elsif( exists($protein2orf{$key}) ) {
                       		print FH "$protein2orf{$key}\t$match_ratio_result{$key}\n";
			}
			else {
                       		print FH "$key\t$match_ratio_result{$key}\n";
			}
		}

		close FH;
		
	}

}


###################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     compute_interpro_similarity 
# PURPOSE    compute_interpro_similarity between two proteins
# ARGS       $protein1 : protein index 1 => this is the base protein
# 	     $protein2 : protein index 2 
# RETURN     $match_ratio 
###################################################################################
sub compute_interpro_similarity {

	my ($protein1, $protein2) = (shift, shift);
	my ($entry_no, $match_no) = (0, 0);

	for(my $i=0; $i<$interpro_cnt; $i++) {
		# currently we check the match of terms in the base protein
		if($protein_interpro_2d_array->[$protein1][$i]==1) {
			$entry_no++;
			
			if($protein_interpro_2d_array->[$protein2][$i]==1) {
				$match_no++;
			}
		}

	}

	if($entry_no>0) {
		return ($match_no/$entry_no);
	}
	else {
		return 0;
	}

}

###################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     init 
# PURPOSE    init the global variables for parsing ProteinAC IDA from file
# ARGS       
# RETURN     etc. initialize the global variable 
#                 @protein_interpro_2d_array,
#                 @protein_vector  
###################################################################################
sub init {
        print "\n\n********************************************\n";
        print "init() ....\n";
        print "********************************************\n";

        #initialize the directory
        &init_dir();

	#initialize the interpro vector from the IDA term file
        print "parse_IDA_term_file() ....\n";
	&parse_IDA_term_file(IDA_TERM_FILE_NAME);

	#initialize the protein vector
        print "init_protein_vector() ....\n";
	&init_protein_vector(PROTEIN_IDA_FILE_NAME);

	#count the number of proteins and the interpro entry 
	$protein_cnt = $#protein_vector + 1; 
	$interpro_cnt = $#interpro_vector + 1;

	#initialize the protein_interpro_2d_array
	$protein_interpro_2d_array = [];
	for(my $i=0; $i<$protein_cnt; $i++) {
		for(my $j=0; $j<$interpro_cnt; $j++) {
			$protein_interpro_2d_array->[$i][$j] = 0;
		}
	}

	# parse the protein_vector, protein_interpro_array  into 2-dimensional array 
        print "parse_ProteinAC_IDA_term() ....\n";
	for(my $i=0; $i<$protein_cnt; $i++) {
		&parse_ProteinAC_IDA_term($protein_interpro_array[$i], $i);
	}	

	# store that result into file 
	if(MY_DEBUG) {

	open FH, ">" . VECTOR_FILE_NAME
		or die "[ERROR] can't open the file " . VECTOR_FILE_NAME . " : $!";
	
	for(my $i=0; $i<$protein_cnt; $i++) {

		print FH "$protein_vector[$i]\t";

		for(my $j=0; $j<$interpro_cnt; $j++) {
			print FH "$protein_interpro_2d_array->[$i][$j] ";
		}

		print FH "\n";
	}

	close FH;

	}

	# initialize the protein2acf hash
	&load_protein2orf(PROTEIN2YEAST_FILE_NAME); 

	# check
	if(MY_DEBUG) {
		for(my $i=0; $i<$protein_cnt; $i++) {
			print "$protein_vector[$i] : " ;
			for(my $j=0; $j<$interpro_cnt; $j++) {
				if($protein_interpro_2d_array->[$i][$j]==1) {
					print "$interpro_vector[$j]\t" ;
				}
			}
			print "\n";
		}
	}
}

###################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     init_dir 
# PURPOSE    make the directory if not exist 
###################################################################################
sub init_dir() {

    if(-e RESULT_DIR) {
	# fine 
    }
    else {
        system "mkdir -p " . RESULT_DIR; 
        print "\nDirectory " . RESULT_DIR . " is generated!\n";
    } 

}

###################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     init_protein_vector 
# PURPOSE    init @protein_vector  from file
# ARGS       $filepath 
# RETURN     etc. initialize the global variable 
#                 @protein_vector  
###################################################################################
sub init_protein_vector {

	my $filepath = shift;
	my ($index, $line, @protein);

	# open the file
	open FH, "<$filepath"
		or die "[ERROR] can't open the file $filepath : $!";

	# read the content
 	$index = 0;	
        while( defined($line=<FH>) ) {
                chomp($line);

                # if this is the comment or blank line, skip it.
		# 5 is determined arbitrary !
                if( $line =~ /^#/  || length($line)<5 ) {  
			next;
		}

		# this is the ProteinAC \t IDA term list
		@protein = split /\t/, $line;

		# use only the ProteinAC at this time
		$protein_vector[$index] = $protein[0];         # ProteinAC
		$protein_interpro_array[$index] = $protein[1]; # IDA List	

		$index++;	
        } 
                                                                                                                             
        # close the file
        close FH;


	if(MY_DEBUG) {
		for(my $i=0; $i<$index; $i++) {
			print "$protein_vector[$i]=>$protein_interpro_array[$i]\n";
		}
		print "\nwhole $index protein interpro pairs are stored\n";
	}

}


#########################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     parse_ProteinAC_IDA_term and register it into the global variable %interpro_has
# PURPOSE    parse ProteinAC IDA terms from file
# ARGS       $termlist   ex. "IDA1025,11011.1965,1005,11011.1965[1841],11011.1965[1841]";
#            $protein_index     
# RETURN     etc. all vector terms are stored the global variable #interpro_vector 
#########################################################################################
sub parse_ProteinAC_IDA_term {

	my ($termlist, $protein_index, $output_path) = (shift, shift, shift);
	my (@terms, $term, @cterms, $cterm, $index);

	#print "$termlist\n";	

	#remove the header IDA
	$termlist = substr($termlist, 3, length($termlist)-3);

	#split the termlist as terms
	@terms = split /,/, $termlist;

	#split the concatenated term
  	foreach $term (@terms) {
		# check if the term is composite or not
		$index = index($term, ".");
		if($index>0) {
			# this is the composite term, split it.
			@cterms = split /./, $term;

			# for every composite elements,
			# register it.
			foreach $cterm (@cterms) {
				$cterm = &clean_IDA_term($cterm);
				&set_protein_interpro_vector($protein_index, $cterm);
			}

		}
		else {
			# this is a simple term. 
			# if this term is not registered yet, 
			# register it with only Interpro Entry.
			$term = &clean_IDA_term($term);
			&set_protein_interpro_vector($protein_index, $term);
		}
	} 	


}

###################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     set_protein_interpro_vector 
# PURPOSE    if a protein has the same term as ARG, set the vector 1 
# ARGS       $protein_index
#            $interpro_term
###################################################################################
sub set_protein_interpro_vector {

	my ($protein_index, $term) = (shift, shift);

	# found the index of $term in the @interpro_vector and 
	# set the correspondencing 2d-array to 1
	for(my $i=0; $i<$interpro_cnt; $i++) {
		# print "[$interpro_vector[$i]] eq [$term]\n";

		if($interpro_vector[$i] eq $term) {
			$protein_interpro_2d_array->[$protein_index][$i] = 1;
			return;
		}
	}

	# something is wrong
	die "[ERROR] no vector term like $term\n";
}

###################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     parse_IDA_term_file
# PURPOSE    parse IDA terms from file
# ARGS       $filepath   ex. "/home/my_perl/IGA/interpro_data/IDA_term.txt";
# RETURN     etc. all vector terms are stored the global variable 
#                 @interpro_vector 
###################################################################################
sub parse_IDA_term_file {

	my $filepath = shift;
	my ($line, @interpro_hash);

	# open the file
	open FH, "<$filepath"
		or die "[ERROR] can't open the file $filepath : $!";

	# read the content
        while( defined($line=<FH>) ) {
                                                                                                                             
                chomp($line);
                                                                                                                             
                # if this is the comment or blank line, skip it.
		# 5 is determined arbitrary !
                if( $line =~ /^#/  || length($line)<5 ) {  
			next;
		}

		# this is the IDA term list
		if ( $line =~ /^IDA/ ) {
			# remove the IDA header
                	&parse_IDA_term( substr($line, 3, length($line)-3) );                                                     
		}
		else {
			warn "[WARNING] strange line $line : $!";
		}

        } 
                                                                                                                             
        # close the file
        close FH;

      	# sort the summary according to the key decreasing order
       	@interpro_vector = sort sort_by_key keys %interpro_hash;

	if(MY_DEBUG) {
		my ($key, $value, $no);
		$no=0;
                foreach $key (@interpro_vector) {
                        print "\t[$no] $key\n";
			$no++;
		}
	}	
} 


#########################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     parse_IDA_term and register it into the global variable %interpro_has
# PURPOSE    parse IDA terms from file
# ARGS       $termlist   ex. "IDA1025,11011.1965,1005,11011.1965[1841],11011.1965[1841]";
# RETURN     etc. all vector terms are stored the global variable #interpro_vector 
#########################################################################################
sub parse_IDA_term {

	my $termlist = shift;
	my (@terms, $term, @cterms, $cterm, $index);

	#print "$termlist\n";	

	#split the termlist as terms
	@terms = split /,/, $termlist;

	#split the concatenated term
  	foreach $term (@terms) {
		# check if the term is composite or not
		$index = index($term, ".");
		if($index>0) {
			# this is the composite term, split it.
			@cterms = split /./, $term;

			# for every composite elements,
			# register it.
			foreach $cterm (@cterms) {
				$cterm = &clean_IDA_term($cterm);
				$interpro_hash{$cterm} = 0 if !exists($interpro_hash{$cterm});
			}

		}
		else {
			# this is a simple term. 
			# if this term is not registered yet, 
			# register it with only Interpro Entry.
			$term = &clean_IDA_term($term);
			$interpro_hash{$term} = 0 if !exists($interpro_hash{$term}) ;
		}
	} 	


}

###################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     clean_IDA_term
# PURPOSE    remove currently unnecessary trailing information for IDA entry 
# ARGS       $term   ex. "1965x^[1841]"
# RETURN     $clean_term   ex. "1965" 
###################################################################################
sub clean_IDA_term {

	my $term = shift;
	my $clean_term;

	my ($index1, $index2);

	$index1 = index($term, "x");
	$index2 = index($term, "[");

	if($index1==-1 && $index2==-1) {
		# this one is clean
		$clean_term = $term;
	}
	elsif($index1==-1) {
		$clean_term = substr($term, 0, $index2);
	}
	elsif($index2==-1) {
		$clean_term = substr($term, 0, $index1);
	}	
	elsif($index1<$index2) {
		$clean_term = substr($term, 0, $index1);
	}
	else {
		$clean_term = substr($term, 0, $index2);
	}

	# debug
	# print "$term=>$clean_term\n";

	return $clean_term;
}

sub sort_by_key {
	$a cmp $b
}

sub sort_by_match_ratio {
	$match_ratio_result{$a} cmp $match_ratio_result{$b};	
}
 
sub init_match_ratio_result {
	my ($protein, $match_ratio);	

	while( ($protein, $match_ratio) = each %match_ratio_result ) {
		$match_ratio_result{$protein} = -1; 
	}
}
